/*
** Modified by Neal May/June 2011
** Remodified by Neal Sept 2014 for Letters piece
** Stripped out lots of things to just compare Logit and clogit for standard normal x
** and fixed earlier version to compute RMS error and bias
** This program is a modified version of
** 
** Bias in Conditional and Unconditional Maximum Likelihood Estimation
**  Ethan Katz (ekatz@fas.harvard.edu)
**  last modified 02/06/01
** of Tom Coupe's
** date: 02 02 2005
** this modification solves the error in Katz's original program as it does allow to include fixed effects
** Tom Coup?  tcoupe@eerc.kiev.ua
** NB - 2014 The modification is now almost 100% so I am no longer calling this a modification 
** NB 2015 This particular code produces Table 2
*/

version 14
clear all
/*This is done once and for all including setting any parms used and controlling looping*/
set matsize 2000
set seed 12947
set more off
capture log close
log using biasnofe2.log, replace
/*the stuff below may change for different sets of runs*/
   matrix times = (5\7\10\20\30\50\75\100\3)
local ntimes=rowsof(times)
matrix groups = (52\50\100)
local ngroups=rowsof(groups)
matrix conterm = (0)
  local nnovar=0 /*number of units where all y are 0*/
local  sims=	1000	/* number of simulations */
 matrix xco=(1)
local nx = rowsof(xco)
local ncons=rowsof(conterm)
local rowres=`ncons'*`ntimes'*`ngroups'*`nx'
matrix define resnew=J(`rowres',7,0)

matrix colnames resnew = N T RMSLOGIT BIASLOGIT  NORESULTPC RMSLOGITJUSTX BIASLOGITJUSTX
noisily display    " SIMS " `sims'
/*nothing changed below*/
local counter=0
  /* Now we loop over constants and T */

  forvalues xfn = 1/`nx' {
  forvalues groupnum = 1/`ngroups' {
forvalues connum=1/`ncons' {
forvalues  timenum = 1/`ntimes'  {
quietly {
local counter=`counter'+1
drop _all
set obs 10000
local xcoeff = xco[`xfn',1]
local  time=times[`timenum',1]
 local n=groups[`groupnum',1]
local concoef= conterm[`connum',1] /* coefficient on constant term */
noisily display "sim number  " `counter' " n  " `n'  " t " `time'
local t=`time'			/* number of time periods */
local nomit=`nnovar'*`t'

local nt=`n'*`t'		
local n1=`n'+1
local n2=`n'+2
local t1=`t'+1


/*now we are going to generate the period and group variable*/
  
gen group = ceil(_n/`t') if _n<=`nt'
gen time = mod(_n-1,`t')+1 if _n<=`nt'
gen rangep = . /*80% range of probs*/
gen avt = . /*storing t ratios of clogit */
gen meanskew=. /*skewness in x*/
gen meankurt=. /*kurtosis in x*/
gen meanrsq = . /*rsq of x on fe1*/
gen errlogit = . /*absolute error in b for X coef of logit*/
    gen errlogittrue = . /*absolute error in b for X coef of logit w true dummy*/
gen errlogits = . /* absolute error in b for x coef of clogit*/
    gen coeflogit = . /*b for X coef of logit*/
        gen coeflogits = . /*  b for x coef of clogit*/
            gen errlogitnofe=.
gen coeflogitnofe=.
    gen meanp=. /*mean of the probs*/
      gen groupdrop = . /* groups dropped by clogit */
          gen perrlogit=.
gen perrlogitc=.
gen bad=.
    
local i=1
while `i'<=`sims' {
/* now we are going to draw the true effects which do not vary by group */

    /*we are going to give x some kurtosis by using a mild lognormal */
        replace bad=0 in `i'
    generate x= rnormal() if _n<=`nt'
    generate y=`concoef' + `xcoeff'*x + rnormal(0) if _n<`nt'
forvalues jj=1/50 {
     local vname="v"+strofreal(`jj')
     gen `vname'= rnormal() if _n<=`nt'
       }
/* draws Y from a Beroulli distribution and estimates parameters */

	



                    
capture reg  y x v1-v50  in 1/`nt'
                       
                        if _rc!=0 {
                           
                   replace errlogit=. in `i'
                            replace coeflogit= . in `i'
                            replace bad=1 in `i'
                            capture drop   y   x v*
                            
                       
                            }
else {
        replace errlogit=(_b[x]-`xcoeff')^2 in `i'
replace coeflogit=_b[x] in `i'
  
    }

capture reg y x in 1/`nt'

       if _rc!=0 {
                           
                   replace errlogits =. in `i'
                            replace coeflogits= . in `i'
                            capture drop    y      x v*
                            
                       
                            }
else {
        replace errlogits=(_b[x]-`xcoeff')^2 in `i'
replace coeflogits=_b[x] in `i'
}
capture drop y      x v*

	local i=`i'+1

      }
           
/*this is the end of a single sim */



/*done after all the sims for a given T and concoef are done */

matrix resnew[`counter',1]=`n'  /*number of groups*/


matrix resnew[`counter',2]=`time' /*number of time periods*/

 
 


summarize errlogit, detail
matrix resnew[`counter',3]=sqrt(r(p50))
summarize coeflogit, detail
matrix resnew[`counter', 4] = r(p50)-`xcoeff'
summar bad
matrix resnew[`counter',5]=r(mean)
summarize errlogits, detail
matrix resnew[`counter',6]=sqrt(r(p50))
summarize coeflogits, detail
matrix resnew[`counter', 7] = r(p50)-`xcoeff'


}  /*quietly*/
}  /*xcoef*/
} /* groups */
} /* constant */
} /* time */
/*Done once at the end after all looping is done*/
noisily {
 
matrix list  resnew
}

capture erase resbiasnofe2.csv
mat2csv, matrix(resnew) saving(resbiasnofe2)

